IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. ?, NO. ?, 



1 



Implementation strategies for hyperspectral 
unmixing using Bayesian source separation 

Frederic Schmidt, Albrecht Schmidt, Erwan Treguier, 
Mael Guiheneuf, Said Moussaoui and Nicolas Dobigeon, Member, IEEE, 



Abstract — Bayesian Positive Source Separation (BPSS) is a 
useful unsupervised approach for hyperspectral data unmixing, 
where numerical non-negativity of spectra and abundances has 
to be ensured, such in remote sensing. Moreover, it is sensible to 
impose a sum-to-one (full additivity) constraint to the estimated 
source abundances in each pixel. Even though non-negativity 
and full additivity are two necessary properties to get physically 
interpretable results, the use of BPSS algorithms has been 
so far limited by high computation time and large memory 
requirements due to the Markov chain Monte Carlo calculations. 
An implementation strategy which allows one to apply these 
algorithms on a full hyperspectral image, as typical in Earth 
and Planetary Science, is introduced. Effects of pixel selection, 
the impact of such sampling on the relevance of the estimated 
component spectra and abundance maps, as well as on the 
computation times, are discussed. For that purpose, two different 
dataset have been used: a synthetic one and a real hyperspectral 
image from Mars. 

Index Terms — Hyperspectral imaging, source separation, 
Bayesian estimation, implementation strategy, computation time. 



I. Introduction 

In visible and near infrared hyperspectral imaging, each 
image recorded by the sensor is the solar light reflected and 
diffused back from the observed planet surface and atmo- 
sphere at a particular spectral band. Under some assumptions 
related to surface and atmosphere properties - e.g., Lam- 
bertian surface, no intimate mixture, no diffusion terms in 
the atmosphere, homogeneous geometry in the scene - each 
measured spectrum - i.e., each pixel of the observed image 
for several spectral bands - is modeled as a linear mixture 
of the scene component spectra (endmembers) QJ-J3). In this 
model, the weight of each component spectrum is linked to 
its abundance in the surface area which corresponds to the 
underlying pixel. The main goal of hyperspectral unmixing 
is to identify the components of the imaged surface and to 
estimate their respective abundances (4|, (5|. 

By considering P pixels of an hyperspectral image acquired 
in L frequency bands, the observed spectra are gathered 
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in a P x L data matrix X, potentially ignoring spatiality. 
Each row of this matrix contains a measured spectrum at 
a pixel with spatial index p = 1, . . . , P. According to the 
linear mixing model, the pth spectrum, 1 < p < P, can 
be expressed as a linear combination of R pure spectra of 
the surface components. Using matrix notations, this linear 
spectral mixing model can be written as 

X^AS (1) 

where non-negative matrices A G if H and S G if xi 
approximate X G lf xi in the sense that |||AS - X|| 2 
is minimized (W*' denotes the space of matrices with only 
non-negative entries of respective dimensions). The rows of 
matrix S now contain the pure surface spectra of the R 
components, and each element a pr of matrix A corresponds to 
the abundance of the rth component in pixel with spatial index 
p. For a qualitative and quantitative description of the observed 
scene composition, the estimation problem consists of finding 
matrices S and A that allow one to explain the data matrix 
X and have a coherent physical interpretation. This approach 
casts the hyperspectral unmixing as a source separation prob- 
lem under a linear instantaneous mixing model (6). Source 
separation is a statistical multivariate data processing problem 
whose aim is to recover unknown signals (called sources) from 
noisy and mixed observations of these sources (7), (8). 

This problem has been studied in-depth in recent years, 
starting with pioneer work more than 15 years ago (9j, |T0| . 
From a statistical point of view, the problem is also related to 
Principal Component Analysis (PCA) and k-means clustering 
(see (TTJ for an overview). Also note that the factorization AS 
is not uniquely defined. For instance, for any matrices Z G 
R% x - such that ZZ 1 = I, then AZZ^S = (AZ)(Z- X S) = 
A'S' is a solution as well; this holds even if the minimization 
is able to find a global minimum. However, when solving this 
separation problem with hyperspectral data, several constraints 
can be considered to reduce the set of admissible solutions. 
A first hard constraint is the non-negativity of the elements of 
both matrices S and A since they correspond to pure spectra 
and abundances of the surface components, respectively. A 
second constraint that may be imposed is the sum-to-one 
(additivity) constraint of the abundances. Indeed the abundance 
weights correspond to proportions and therefore should sum 
to unity. 

Several algorithms have been proposed in the literature to 
solve fully constrained unmixing problems, i.e., handling both 
of the constraints imposed on the spectra and abundances. 
Specifically, an iterative algorithm called ICE (Iterated Con- 
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strained Endmembers) has been proposed in (12) to minimize 
the size of the simplex formed by the estimated endmembers. 
However, as noted in (13) , results provided by ICE strongly 
depend on the choice of the algorithm parameters. More re- 
cently, Jia and Qian have developed in (14) complexity-based 
BSS algorithms that exploit the pixel correlations to recover 
endmember signatures. In (T5) , Miao et al. have introduced a 
non-negative matrix factorization (NMF) algorithm with an 
additivity penalty on the abundance coefficients. Similarly, 
other constrained NMF approaches exploiting smoothness and 
sparseness features have been considered in (16). Note that all 
the strategies described above are based on an optimization 
scheme to minimize a penalty criterion. Consequently, they 
may suffer from convergence issues, e.g., due to the presence 
of local maxima and the large number of parameters to be 
estimated. 

Alternatively, the constrained separation problem can be 
conveniently addressed in a Bayesian framework. Two algo- 
rithms that perform unsupervised separation under positivity 
and sum- to-one constraints have been recently proposed (17), 
|T8) . These algorithms are based on hierarchical Bayesian 
modeling to encode prior information regarding the obser- 
vation process, the parameters of interest and include the 
positivity and full additivity constraints. The complexity of the 
inference from the posterior distribution of the parameters of 
interest is tackled using Markov chain Monte Carlo (MCMC) 
methods (19), (20), which has been proposed to analyze 
hyperspectral images (21). The algorithmic details are not 
described here. The reader is invited to read the pseudo- 
codes summarized in Algo.'s 1 and 2, and to consult (T7) 
and (T8) for further details. The only difference between 
the two BPSS algorithms is the sampling of the abundance 
vectors a. p (p = 1, . . . , P). However, since these algorithms 
rely on MCMC methods, the computation time drastically 
increases with the image size and these algorithms have not 
been applied for large scale data processing in spite of their 
high effectiveness. 

The aim of this article is to discuss some implementation 
strategies which allow one to apply these algorithms to real 
hyperspectral data, even if images are large. Previous works 
about blind source separation of hyperspectral images have 
been proposed |22)-|24) but only few use positivity /sum- 
to-unity constraints |25| . To overcome this difficulty, a first 
approach has been proposed in (25) to combine Independent 
Component Analysis (ICA) and Bayesian Positive Source 
Separation (BPSS). Firstly, applying an ICA algorithm (such 
as JADE (26) or FastlCApT)) to the hyperspectral images is 
applied to get a rough spatial classification of the scene and to 
sample relevant pixels (i.e., from each class, the pixels whose 
spectra are mostly uncorrected are selected). Secondly, the 
spectra associated to these pixels will serve in the Bayesian 
separation algorithm to estimate the endmember spectra. Fi- 
nally, the abundances can then be estimated on the whole 
image using the estimated spectra. However, this strategy 
presents a limitation related to the difficulty to determine the 
number of pixels to retain from each independent component 
class. In this paper, another pixel selection strategy based on 
the computation of the convex hull of the hyperspectral data 



is introduced. Its influence on the separation performances is 
also discussed. The issue of estimating the number of sources, 
or "intrinsic dimension" (28), will not be addressed in this 
article. Several methods have been proposed in the literature 
(29), (30). 

This paper is organized as follows. Section [IT| describes the 
proposed implementation strategies adopted for this work. Sec- 
tion [In] summarizes the improvements related to the technical 



aspects of memory storage and computation issues. Section IV 
discusses the performances of the resulting algorithms when 
the pixel selection preprocessing step is introduced. 

II. Optimization Strategies 

The optimization consists of two independent parts which 
will be referred to: (i) Technical Optimization (TO) to reduce 
the memory footprint, the average cost of algorithmic oper- 
ations, and make smarter reuse of memory (ii) Convex Hull 
Optimization (CHO) to reduce the number of spectra to be 
processed. 

Both parts enabled us to analyze hyperspectral images that 
so far were not open to analysis. The authors stress that the 
techniques applied in (i) do not alter the results of the original 
algorithm (see section [III]). On the other hand, the optimization 
strategy (ii) only selects a subset of the original input and 
therefore may change the results. Impact of the strategy (ii) 
needs to be evaluated, which will be presented in section [TV 



A. Technical Optimization (TO) 

The algorithms introduced in (17) , (T8) and referred to 
as BPSS and BPSS2, respectively, could be successfully 
launched on an image of a restricted size, typically of a few 
thousand pixels. The main goal of this work is to optimize the 
memory requirement of these algorithms to process a whole 
hyperspectral image of 100000 spectra as it typically occurs 
in Earth and Planetary Science. Since the time requirements 
of the computation increase drastically for a larger number 
of pixels and a larger number of sources, another challenging 
objective is to reduce as much as possible the computation 
time. In that respect, our proposal is to discuss the memory 
storage, the data representation, the operating system archi- 
tecture and the computing parallelization. These algorithms 
have been implemented in MATLAB© for this work but future 
implementations will be done in other languages as well. 

1) Memory: Thanks to the MATLAB© profiler, it can be 
noticed that the main limitation of the BPSS implementation 
is the contiguous memory. Fragmentation may occur when 
variables are resized after memory allocation. In this case, the 
memory management might not be able to allocate a chunk 
of memory that is large enough to hold the new variable. 
Significant garbage collection may set in, which may have 
a significant performance impact. In our case, to reduce the 
impact of garbage collection, pre-allocating the matrices and 
work with global variables has been found to be useful. 

2) Precision: MATLAB© by default computes on double 
precision. However computing with single data type saves a 
lot of computation time while providing sufficient arithmetic 
precision. It has been estimated to win up to 60% computation 
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Algorithm 1 Bayesian positive source separation algorithm 
(BPSS) 

for i = 1 , . . . , -/Vmc do 

% sampling the abundance hyperparameters 
for p = 1 , . . . , P do 
Draw A p from the pdf 



/ (Ap|a p: ,7 P ) oc Y[ 



r(A p ) ^ 



Ap 1r+(A p ). 



end for 

% sampling the abundance hyperparameters 
for p = 1 , . . . , P do 

Draw 7 P from the gamma distribution 



7 P | A p , a p: ~ £ I 1 + i?Ap + e, ^ a p , r + . 



end for 

% sampling the abundance vectors 
for p = 1 , . . . , P and r — 1 , . . . , R do 
Draw a PjT . from the pdf 

/ (a P;r |Ap, 7p, S, <Tg, X) 



00 a p,r ll R+( a p,r)exp 



end for 

% sampling the noise hyperparameter 

Draw ip e from the inverse-gamma distribution 



tfte CT e , p e ~ X£ 



f P Pe 1 ^ 1 

1 2 ' 2 ^-J a e % 

\ p=i e ;P 



% sampling the noise variances 
for p = 1 , . . . , P do 

Draw <jg p from the inverse-gamma distribution 



cr e v \ip e , a p: , S, x p: ~ X£ 



p e + £ + ||Xp: - Sa p: | 



end for 

% sampling the source hyperparameters 
for r = 1, . . . , R do 
Draw a r from the pdf 



end for 

% sampling the source hyperparameters 
for r = 1, . . . , R do 

Draw /3 r from the gamma distribution 



Algorithm 2 Fully constrained Bayesian positive source sep- 
aration algorithm (BPSS 2) 

for i = 1 , . . . , TVmc do 

% sampling the abundance vectors 
for p = 1 , . . . , P do 
Draw a p: from the pdf 

/ (a p: |A,^,x) 

00 6XP [~2~ ~ M p) T A p 1 ( ap: ~ 1§ ( a P ; )- 



with 



a p: ; a P;r > 0, Vr = 1, . . . , R, a P;r = 1 



end for 

% sampling the noise hyperparameter 

Draw from the inverse-gamma distribution 



I z Pp e 1 ^ 1 

p=l e,p 



% sampling the noise variances 
for p = 1 , . . . , P do 

Draw cr^ from the inverse-gamma distribution 



CT e |-0e,ap:,S,Xp : ~ XQ 



Pe + L 1p e + ||Xp: - Sa p : 



end for 

% sampling the source hyperparameters 
for r = 1 , . . . , R do 
Draw a r from the pdf 



/(a P |. P! ,ft.) ocJl I^T^jl e—-l, + (a r ) 



end for 

% sampling the source hyperparameters 
for r = 1 , . . . , R do 

Draw (3 r from the gamma distribution 



end for 

% sampling the source spectrum 
for r = 1 , . . . , P and I = 1 , . . . , L do 
Draw s r ,z from the pdf 

/ ( s r ,z lav, /3 r , A, <Tg, X) 



^ S r,7 ( S r*,l) ex P 



252 



end for 
end for 



/3 r |a r , s r: ~ Q ( 1 + Lav + e, ^ s r , z + . 



end for 

% sampling the source spectrum 
for r = 1 , . . . , R and I = 1 , . . . , L do 
Draw s r ^ from the pdf 

/ (s r , z |a r ,/3 r ,A,o-;f,x) 



00 S rJ ll M+(«r,z) ex P 



Qr,Z - p r , Z ) 2 
2(52 



■ (3 r s r ,i 



end for 
end for 



time on an x86 processor architecture, while the changes to 
the code have been minimal. Furthermore, most dataset come 
as single precision. 

3) OS Architecture: It is interesting to note that 
MATLAB® is limited in terms of memory usage (regardless of 
the size of physical memory). This depends on the Operating 
System (OS) and on the MATLAB version (see Table |II- A3 } . 

Therefore, a 32-bits LINUX architecture has been chosen. 

4) Parallelization: MATLAB® contains libraries dedicated 
to automatically parallelize parts of the algorithms on a single 
computer. BPSS has been run on a 4-core machine. The 
underlying matrix libraries already provide a certain level 
of parallelism depending on the number of available cores. 
However, in the future, parts of the code could be parallelized 
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lVTpmrvrv T imitation 


32-bit Microsoft Windows XP 
Windows Vista 


2GB 


32-bit Windows XP with 3 GB boot.ini switch 
32-bit Windows Vista with increaseuserva set 


3GB 


32-bit LINUX 


3GB 


64-bit Windows XP, Linux, 
Apple Macintosh OS X 
or SunSolaris running 32-bit MATLAB© 


4GB 


64-bit Windows XP, Windows Vista, 
Linux, or Solaris running 64-bit MATLAB© 


8GB 



TABLE I 

Summary of memory limitation depending on Operating System. 



and the jobs could be submitted to a grid in order to speed up 
the calculation process. 



B. Convex Hull Optimization (CHO) 

The proposed pixel selection strategy is based on the convex 
hull of the data matrix projection into the subspace spanned 
by the principal components. The convex hull of a point set 
is the smallest convex set that includes all the points (3T| . 
The pixels associated to the vertices of the convex hull are 
selected (32) and are expected, despite their limited number, 
to exhibit the main spectral features of the whole dataset. In 
terms of abundances, this sample of points should contain the 
pixels with the highest abundances of the components which 
contribute to the investigated hyperspectral image (i.e., the 
purest pixels or most extreme pixels). It can be used as a 
concise representation of the dataset which still features the 
strongest spectral signatures available in the original image. 
This strategy is also used as a first step in endmember 
extraction algorithms for dimension reduction and purest pixel 
determination p2|-|37|. Pixel selection has the advantage, 
to reduce the number of mixture spectra to unmix and to 
enforce the sparsity of the mixing coefficients to be estimated. 
Note that the spectral dimension of the selected spectra is not 
changed, only the spatial dimension is reduced since only few 
pixels are selected. 

The convex hull selection has been implemented after 
seven spectral components have been selected through PCA, 
which turned out to be a good compromise between resource 
consumption and accuracy. 

III. Performance and Accuracy of TO 

All the following runs are performed on a Quad-Core AMD 
Opteron(tm) Processor 8384 at 2.7GHz with 2Gb of memory. 

1 ) Performance: Computation times between the previous 
version of BPSS and the TO version have been compared when 
processing a synthetic dataset of 1052 spectra of 128 bands 
and 3 sources. For a run attempting at estimating 3 sources, the 
computation time has decreased from 1106s (previous version) 
to 724s (TO version), i.e., by a factor of about 1.5. In addition, 
the total memory consumption is nearly half for the TO version 
of the algorithm. 



2) Accuracy: Due to the stochastic nature of the BPSS 
algorithms, it is difficult to demonstrate that two algorithms 
are semantically identical. In order to check that no significant 
loss of accuracy has been induced by the TO and especially 
by the change from double to single precision, several tests 
have been performed with different random seeds xi an d X2> 
used for the initialization step of the MCMC. The sources S x i 
and S X 2 estimated with and without TO have been compared. 

The average correlations between S x i and S X 2 without TO 
are 0.9816 ± 0.0315and 0.9818 ± 0.0255 with TO. These 
correlations are due to the stochastic approach in the Bayesian 
framework. Correlation values are similar, indicating that the 
stochastic variance has not been affected by TO. 

The average cross-correlation between S x i and S X 2 with 
and without TO is 0.9760 ± 0.0388. This value is similar to 
the correlation due to stochastic process, demonstrating that 
the TO version is equivalent to the original version of BPSS. 

No significant differences have been observed, confirming 
that the TO version is equivalent to the original version of 
BPSS. 

IV. Performance and Accuracy of CHO 

The impact of the convex hull pixel selection pre-processing 
step has been evaluated on two dataset: (i) synthetic data 
generated from linear mixtures of known materials and (ii) an 
OMEGA hyperspectral image of the south polar cap of Mars 
as an example from Planetology. Since the BPSS with TO 
has been shown to be semantically equivalent to the previous 
version, the TO approach is used in the rest of this article. 

A. Synthetic data 

1) Description: Several synthetic dataset have been gen- 
erated by mixing a known number of endmembers, with 
abundances simulated with uniform distribution. The generated 
dataset are of size 200 x 500 pixels, which is a spatial size 
similar to the one of a typical hyperspectral image. For the 
endmembers, the following spectra have been used: H 2 and 
CO2 ice spectra (38), (39) and mineral spectra from the USGS 
Digital Spectral Library splib06a (40), resampled to match 
the 128 wavelengths of OMEGA C Channel j?T). To ensure 
the sum- to-one constraint on the R endmember abundances, a 
uniform distribution on the simplex has been used following 
a well established scheme (42). Synthetic dataset have been 
generated with 3, 5 and 10 endmembers. Based on this method, 
dataset for which the maximum abundance of each single 
endmember was limited to a certain value (100%, 80% and 
60%) have also been considered. This latter data, that are 
called "cutoff in the sequel, allows one to test the method 
efficiency face to various conditions in terms of purity of 
the samples (in cases where pure - to a certain degree - 
components occur in the dataset or not). In addition, a 3 
component asymmetric dataset has been investigated, with one 
of the the component abundance (albite) being limited to a 
cutoff of 35% and the abundances of the two others (ices) not 
being limited. Besides, dataset with some added OMEGA-like 
Gaussian noise, amplified or not, have been also generated and 
investigated. The noise estimation on the dark currents of the 
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Algorithm 


Without CHO 


With CHO 


Time ratio 


BPSS 


71463 (BPSS-2) 


2205 (BPSS-1) 


32.41 


BPSS2 


530654 (BPSS2-2) 


7133 (BPSS2-1) 


74.39 



TABLE II 

Computation times (after TO) in seconds, for a synthetic 

DATASET WITH 3 ENDMEMBERS (NO CUTOFF, NO NOISE), FOR BOTH BPSS 
AND BPSS2, WITH AND WITHOUT CHO. IN THIS EXAMPLE, 944 PIXELS 
WERE SELECTED FOR THE CHO, AM ONG A TOTAL OF 100000. THE NAME 
OF THE RUN OF TABLES |lV- A3 G| AND llVllS NOTED IN PARENTHESIS. 



OMEGA instruments for observation 41_1 has been used 139) . 
Note that for all the considered simulation scenarii, the number 
of sources to be estimated has been tuned to the actual number 
of endmembers used to produce the artificial dataset. 

2) Performance: Computation times are about 50 times 
shorter when pixel selection by convex hull (CHO) is per- 
formed as a preprocessing step (see Table [II]). 

3) Accuracy: 

a) Analysis of the results: The spectrum of each esti- 
mated source has been compared to the spectra from the spec- 
tral library containing the pure endmembers used to produce 
the synthetic dataset. The absolute value of the correlation 
has been used as a similarity measurement, thus as a criterion 
for the determination of the best spectral match. On Figures [T] 
each source is represented along with its best match, according 



to the aforementioned criterion. Table |IV-A3g| (resp. Table [rV) 
shows the results for BPSS (resp. BPSS2). 

A source is considered a good estimation of a certain 
endmember if both are each other best spectral match and 
if their absolute correlation is greater than 80%. For each run, 
the number of well-estimated sources is mentioned in Tables 



IV-A3g| and [IV] Note that endmembers matched by several 
sources, in case it happens, are only counted once. Along 
with the number of well-estimated sources, the mean value 
of the correlations between (only) the well-estimated sources 
and their best spectral match also helps to the assessment of 
the accuracy for the estimation of the whole set of sources for 
each run. Simple distance could not be used here because the 
scale in usual blind source separation is undetermined |8|. 

b) BPSS vs. BPSS2: In most of the tested cases, the 
quality of the estimation is unambiguously better with BPSS2 
than with BPSS (see Tables [TV-A3g| and |IV|). The improvement 



appears to be even more significant when the number of 
endmembers is increasing. Our 3 endmember test dataset is 
a mixture of two endmembers with strong spectral signatures 
(CO2 ice and H2O ice) and a third one with weaker signatures 
(albite), as often with minerals. Interestingly, while using 
BPSS allows one to correctly estimate the ices spectra but 
not albite, BPSS2 is actually able to correctly estimate the 
three endmembers. This confirms that adding the sum-to-one 
constraint is necessary when dealing with such dataset, which 
is important regarding the analysis of other dataset. 

c) Effect of the pixel selection (CHO): With the excep- 
tion of the asymmetric dataset (see below), the endmembers 
is less well-estimated when a pixel selection has been per- 
formed, the loss seeming less significant when the number of 
endmembers is low. 

Also note that the results with pixel selection do not appear 



to be very sensitive to the cutoff variations: the loss of quality 
(between runs performed with and without pixel selection) is 
similar for cutoffs of 60%, 80% and 100%, which can be 
explained by the pixel selection's ability to extract the purest 
available pixels. 

d) Effect of the number of endmembers: Due to curse of 
dimensionality, the more endmembers to be estimated with 
the fixed number of wavelength, the more difficult is the 
estimation |43| , (44). Still, BPSS2 gives excellent results even 
for 10 sources, as all spectra are estimated with a correlation 
coefficient higher than 99% (see fig. [T}. 

e) Effect of the maximum abundance cutoff: The cutoff 
affects the quality of the estimation, which is clearly better, 
for BPSS and BPSS2, when pure components occur in the 
dataset. This has to be remembered when dealing with real 
dataset. 

f) Effect of noise: The results clearly show that the 
method is very robust to noise, as the estimation of the sources 
does not appear to be significantly affected by the addition 
of a Gaussian OMEGA-like noise to the synthetic dataset. 
BPSS2 (without pixel selection) even manages to successfully 
overcome the addition of a 100-times amplified OMEGA-like 
noise (see Tables [TV] and [2]). 

g) Effect of asymmetry in maximum abundance cutoff : 
In this case, the results are better with pixel selection rather 
than without. BPSS 2 with pixel selection is the only run 
(performed on this synthetic dataset) that allows one to suc- 
cessfully estimate the three endmembers that have been used to 
generate the dataset, including albite, whose abundances have 
been limited to a cutoff of 35% and whose spectral signature 
is weaker than the ones of the other endmembers (ices). This 
result can be explained by the fact that pixel selection is able to 
extract the pixels with the strongest available albite signature, 
and consequently overcomes the blinding effect of the ices 
occurring in the whole dataset, that has affected the results 
when no pixel selection has been performed. 




Fig. 1. Sources estimated by BPSS2 (blue lines) and their spectral matches 
(red dotted lines), for an artificial dataset with 10 endmembers (no cutoff, no 
noise). 
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Run Id 


Cutoff (%) 


Nb of endmembers 


Noise 


Pixel selection 


Nb of well-estimated sources 


Mean correlation (%) 


BPSS-1 


100% 


3 


no 


yes 


2/3 


on n /i a 1 
99.9441 


one c o 


100% 


3 


no 


no 


2/3 


99.9963 




80% 


3 


no 


yes 


2/3 


no coo/; 
95.5826 


BPSS-4 


80% 


3 


no 


no 


2/3 


99.9286 


BPSS-5 


60% 


3 


no 


yes 


2/3 


no one 

98.2135 


rSroo-O 


OUvo 


J 


no 


no 


LI D 


QQ O/IQQ 


nnc c h 


100% 


5 


no 


yes 


2/5 


92.0689 


RPQC O 

oroo-o 


IUUto 


c 

J 


no 


no 




no ci oi 
yL.j lol 




100% 


10 


no 


yes 


5/10 


o o i no 1 
88.3U81 


T3DQQ IH 


IUUto 


IU 


no 


no 


O/IU 


Q/L 


T3DQQ 1 1 


IUUto 




U1V1 11 Li A 


yes 


LI D 


QQ QCifTl 

yy.yyj'o 1 


T3DQQ 10 


IUUto 


-2 


U1V1 HO A 


no 


LI D 


QQ QQAT 


one c ii 


IUUto 




I UXUMLCr A 


yes 


111 
LI D 


QQ QQ07 

99.892 / 


BPSS-14 


100% 


3 


lOxOMEGA 


no 


2/3 


99.9976 


BPSS-15 


100% 


3 


lOOxOMEGA 


yes 


2/3 


98.8182 


BPSS-16 


100% 


3 


lOOxOMEGA 


no 


2/3 


98.2574 


BPSS-17 


ices: 100%, alb.: 35% 


3 


no 


yes 


2/3 


99.9811 


BPSS-18 


ices: 100%, alb.: 35% 


3 


no 


no 


2/3 


98.9855 



TABLE III 

Results obtained for different synthetic dataset with the BPSS algorithm. Characteristics of each dataset are shown: number 

OF ENDMEMBERS, CUTOFF, AND NOISE. EACH DATASET HAS BEEN ANALYZED WITH A NUMBER OF SOURCES TO BE ESTIMATED EQUAL TO THE NUMBER 
OF ENDMEMBERS USED TO GENERATE THE ARTIFICIAL DATASET, WITH AND WITHOUT PIXEL SELECTION. QUALITY OF THE ESTIMATION IS EXPRESSED 
THROUGH THE NUMBER OF WELL-ESTIMATED SOURCES AND THE MEAN ABSOLUTE EXPRESSION AS EXPLAINED IN THE TEXT. 



Run Id 


Cutoff (%) 


Nb of endmembers 


Noise 


Pixel selection 


Nb of well estimated sources 


Mean correlation (%) 


BPSS2-1 


100% 


3 


no 


yes 


3/3 


99.8923 


BPSS2-2 


100% 


3 


no 


no 


3/3 


99.9222 


BPSS2-3 


80% 


3 


no 


yes 


2/3 


95.8934 


BPSS2-4 


80% 


3 


no 


no 


3/3 


99.9200 


BPSS2-5 


60% 


3 


no 


yes 


2/3 


95.2965 


BPSS2-6 


60% 


3 


no 


no 


3/3 


97.5408 


BPSS2-7 


100% 


5 


no 


yes 


3/5 


99.2821 


BPSS2-8 


100% 


5 


no 


no 


5/5 


99.9174 


BPSS2-9 


100% 


10 


no 


yes 


5/10 


98.9439 


BPSS2-10 


100% 


10 


no 


no 


10/10 


99.9535 


BPSS2-11 


100% 


3 


OMEGA 


yes 


3/3 


99.8726 


BPSS2-12 


100% 


3 


OMEGA 


no 


3/3 


99.9955 


BPSS2-13 


100% 


3 


lOxOMEGA 


yes 


3/3 


99.7298 


BPSS2-14 


100% 


3 


lOxOMEGA 


no 


3/3 


99.9962 


BPSS2-15 


100% 


3 


lOOxOMEGA 


yes 


2/3 


95.4706 


BPSS2-16 


100% 


3 


lOOxOMEGA 


no 


3/3 


98.5647 


BPSS2-17 


ices: 100%, alb.: 35% 


3 


no 


yes 


3/3 


95.9402 


BPSS2-18 


ices: 100%, alb.: 35% 


3 


no 


no 


2/3 


99.7202 



TABLE IV 

Results obtained for different synthetic dataset with the BPSS2 algorithm. Characteristics of each dataset are shown: number 

OF ENDMEMBERS, CUTOFF, AND NOISE. EACH DATASET HAS BEEN ANALYZED WITH A NUMBER OF SOURCES TO BE ESTIMATED EQUAL TO THE NUMBER 
OF ENDMEMBERS USED TO GENERATE THE ARTIFICIAL DATASET, WITH AND WITHOUT PIXEL SELECTION. QUALITY OF THE ESTIMATION IS EXPRESSED 
THROUGH THE NUMBER OF WELL-ESTIMATED SOURCES AND THE MEAN ABSOLUTE EXPRESSION AS EXPLAINED IN THE TEXT. 



B. OMEGA data 

1) Presentation: The OMEGA (Observatoire pour la 
Mineralogie, l'Eau, les Glaces et l'Activite) instrument is 
a spectrometer on board Mars Express (European Space 
Agency), which provides hyperspectral images of the Mars 
surface, with a spatial resolution from 300m to 4km, 96 
channels in the visible range and 256 wavelength channels 
in the near infra-red |45| . In this work, 184 spectral bands 
have been selected according to the best signal to noise ratio. 
Conversely, spectral bands that contain the thermal emission 
have been removed. 

Blind source separation on this dataset has been initiated 
by using the JADE algorithm (46). In particular the image 
41_1 of the permanent south polar region has been used for 



supervised classification approach with WAVANGLET 
unsupervised classification approach 
blind source separation using BPSS 



39), 

and unsupervised 
Since no ground 
truth is available, the results from physical non-linear inversion 
have been considered as a reference [39], (49), (50) . In this 
image, the surface is dominated by dust and some spectra 
contains CO2 and water ices (see fig. [3]). This reference dataset 
for hyperspectral classification is available online Q The Luo et 
al. method introduced in [30] has estimated 2 sources for both 
41_1 and 41_1.CUT images. From previous work using band 
ratio detection [41], physical inversion of the radiative transfer 
(49), (50), supervised classification approach using Wavanglet 



1 http://sites.google.com/site/fredericschmidtplanets/Home/ 



hyperspectral- reference 
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CORRELATION = 0.999289 CORRELATION = 0.997649 CORRELATION = 0.960003 






B0 100 120 



Fig. 2. Sources estimated by BPSS2 (blue lines) and their spectral matches 
(red dotted lines), for an artificial dataset with 3 endmembers and 100-times 
amplified OMEGA-like noise (no cutoff). 



p9| and unsupervised classification (47), 3 endmembers have 
been detected: dust, CO2 and water ice. The number of sources 
has been tuned to 3 in our study. 

The proportion of pixels containing CC^ice and H 2 ice 
on the 41_1 image is estimated to be 16.76% and 21.84%, 
respectively (39). The first 300 lines of the 41_1 image (subset 
named 41_1.CUT) contain all spectra containing ices. For this 
subset, the proportion of pixels with detected CO2 and H2O 
is 48.72% and 63.48%, respecively. 



Algorithm 


Without CHO 


With CHO 


Time ratio 


BPSS 


166 400 (OMEGA-5) 


1 468 (OMEGA-7) 


113.28 


BPSS2 


332 680 (OMEGA-6) 


3 176 (OMEGA-8) 


104.71 



TABLE V 

Computation times in seconds, for OMEGA 41_1 image with 3 

ENDMEMBERS, FOR BOTH BPSS AND BPSS2, WITH AND WITHOUT CHO. 

In this example, 670 pixels have been selected with t he CHO , 

AMONG A TOTAL OF 111488. THE NAME OF THE RUN OF TABLE [rV-B3cl lS 
NOTED IN PARENTHESIS. 



Run Id 


Image 


Algo 


pixel 
selection 


H 2 


C0 2 


dust 


OMEGA- 1 


41_1.CUT 


BPSS 


no 


0.883 


0.955 


0.542 


OMEGA-2 


41_1.CUT 


BPSS2 


no 


0.823 


0.958 


0.980 


OMEGA-3 


41_1.CUT 


BPSS 


yes 


0.956 


0.951 


0.766 


OMEGA-4 


41_1.CUT 


BPSS2 


yes 


0.894 


0.910 


0.975 


OMEGA-5 


41_1 


BPSS 


no 


0.773 


0.957 


0.555 


OMEGA-6 


41_1 


BPSS2 


no 




0.953 


0.512 


OMEGA-7 


41_1 


BPSS 


yes 


0.940 


0.953 


0.372 


OMEGA-8 


41_1 


BPSS2 


yes 


0.450 


0.954 


0.982 



TABLE VI 

Results on algorithms BPSS and BPSS2 on a portion of OMEGA 

IMAGE (41_1.CUT) AND ON THE ENTIRE IMAGE (41_1). FOR 41_1.CUT, 
THE PROPORTION OF PIXELS WITH DETECTED CO2 IS 48.72% AND 
RESPECTIVELY 63.48% FOR H 2 (39). FOR 41_1, THE PROPORTION OF 
PIXELS FOR C0 2 AND H 2 IS 16.76% AND 21.84%. THE COLUMNS H 2 0, 
C0 2 AND DUST INDICATE THE CORRELATION COEFFICIENT BETWEEN THE 
ESTIMATED SOURCES AND THE REFERENCE SPECTRA. (-) INDICATES THAT 
NO IDENTIFICATION OF H 2 NEITHER FROM SPECTRAL NOR SPATIAL 

RESULTS. This source has been detected to BE CO2 ice 

(CORRELATION 0.911). 



Spectral Library Plot 




100 .150 
Band index 



250 



Fig. 3. Reference spectra of the OMEGA hyperspectral image 41_1: (i) in 
blue: synthetic H2O ice with grain size of 100 microns, (ii) in red: synthetic 
CO2 ice with grain size of 10 centimeters, (iii) in black: OMEGA typical 
dust materials with atmosphere absorption. 

2) Performance: Computation times are about 100 times 
shorter when pixel selection by convex hull (CHO) has been 
performed as a preprocessing step (see Table [V). 

3) Accuracy: Table IV-B3c reports the results from differ- 
ent tests, each run is defined by a number. To estimate the 
quality of the estimation, the correlation between the reference 
spectra and the estimated sources has been computed. The 
attribution of each source has been done ad hoc using both 
spectral source and spatial abundances. 

a) Asymmetric abundances of the sources: The quality of 
estimation with both BPSS and BPSS2 is significantly lower 



for dataset 41_1 (run OMEGA-5 to OMEGA-8) in comparison 
with 41_1.CUT (run OMEGA- 1 to OMEGA-4). This result 
suggests that both BPSS and BPSS2 are less efficient in a 
case of an asymmetric distribution of the sources. 

b) BPSS vs. BPSS2: The algorithm BPSS gives signifi- 
cantly better results than BPSS 2 (for instance run OMEGA-3 
vs OMEGA-4). This is due to non-linearity in the radiative 
transfer and noise in the dataset in contradiction with the full 
additivity constraint. 

c) Effect of the pixel selection: When the convex hull se- 
lection has been used as a pre-processing step to BPSS/BPSS2, 
the estimation is significantly better (see fig. |4j for run 
OMEGA-5 and fig. |5] for run OMEGA-7). These results show 
that pixel selection is a way to better take into account the 
occurrence of rare endmembers and thus is an interesting 
method to provide better results. 

V. Discussion and conclusion 

For the first time, a MCMC-based blind source separation 
strategy with positivity and sum-to-one constraints has been 
effectively applied on a complete hyperspectral image with 
a typical size frequently encountered in Earth and Planetary 
Science. The optimization of BPSS (17) and BPSS2 (18) 
presented in this article consists of two independent parts: (i) 
Technical Optimization (TO) reduces the memory footprint, 
lowers the average cost of algorithmic operations, and makes 
smart re-use of memory (ii) Convex Hull Optimization (CHO) 
reduces the number of spectra to process. 
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endmember with abundances until 100%) have demon- 
strated that the sources estimated by the TO strategy 
is equivalent to the CHO strategy (for instance: runs 



Fig. 4. Estimation of 3 sources of the entire OMEGA image 41_1 with BPSS 
using a preprocessing step of pixel selection using the convex hull method. 
The first and third source are clearly identified as CO2 and H2O ices (see 
fig. [3} with a correla tion coefficient of 0.953 and 0.940 (see run OMEGA- 
7 of Table |IV-B3c) . The spatial abundances is well estimated regarding 
the WAVANGLET classification method (39) , [48] . The second source is 
identified as dust with a lower correlation coefficient (0.372). 



ORB0041J 


_V03.REF_run3 


0.25 - 










0.2- 










0.15 - 








/ 0.08 


0.1 / 








/ °-° 6 - 
' 0.04 


0.05- 








0.02 
















Fig. 5. Estimation of 3 sources of the entire OMEGA image 41_1 with 
BPSS without pixel selection. The second source is clearly identified as CO2 
ice (s ee fig. |3) with a correlation coefficient of 0.957 (see run OMEGA-5 of 
Table |IV-B 3c) . The first and third sources are identified to dust and water ice 
with lower correlation coefficients of 0.555 and 0.773. The spatial abundances 
of water ice is not well estimated regarding the WAVANGLET classification 
method (39). 



Figure [6] summarizes the following results in a schematic 
form. 

1) The TO, for both BPSS and BPSS2, allows one to 
decrease the computation times by a factor of 1.5, 
without altering the accuracy of the results. Memory 
consumption has also been reduced by a significant fac- 
tor. With such unambiguous advantages, the TO versions 
of BPSS and BPSS2 can be rather used than the original 
implementations. 

2) Trivially, results obtained for linear artificial dataset 
(with uniform abundance distributions identical for each 



BPSS-1 to BPSS-2 in Table IV-A3g and runs BPSS2-1 



3) 



to BPSS2-2 in Table IV). In this case, pixel selection is 
still relevant to reduce the computation time about 50 
times (Table |TT|). 

Results obtained for artificial dataset with uniform abun- 
dance distributions and identical cutoffs for all endmem- 
bers have shown that the estimation of the sources is 
less accurate when a pixel selection (CHO) has been 



performed (runs BPSS-3 to BPSS-6 in Table |IV-A3g 
and runs BPSS2-3 to BPSS2-6 in Table[TV]). In this case, 
despite 50 times shorter computation times, using pixel 
selection as a preprocessing step seems to be inadequate. 

4) For OMEGA data, the computation time reduction due 
to CHO has been around 100 (Table [V]). Abundance 
distributions can be significantly unbalanced (some end- 
members are significantly less present in the scene). In 
that case, pixel selection by convex hull (CHO) is a 
way to overcome the bias caused by the overwhelming 
endmembers. This has been supported by the results 
obtained for the synthetic dataset of linear mixture using 
unbalanced uniform distribution. 

5) BPSS 2 seems to better estimate the sources in the 
artificial dataset but not in the real case. This is probably 
due to non-linearity or non-Gaussian noise effect. 

6) The method BPSS2 appears to be very robust to Gaus- 
sian noise, as shown by the results obtained on synthetic 
dataset, even with 100 times actual OMEGA noise. 

7) Sometimes, some sources have been well estimated but 
anti-correlated with the real spectra. This behavior has 
been interpreted to be due to linear dependent endmem- 
bers. In that case, spectra built by a linear combination of 
all sources except the considered source already contain 
spectral signatures of the considered source. The last 
source is then anti-correlated with the corresponding 
endmember to decrease his contribution. This behavior 
has to be studied in further details because it is clearly 
a limitation of blind source separation. 

In the future, the choice of the number of sources, which is 
an input in the current implementation, should be automated 
to allow one batch processing without human intervention. A 
methodology of pixel selection for use across dataset should 
also be established to enables integration of source separation 
techniques into larger systems and aim at the generation of 
catalogs and maps. 
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Symmetric abundances 



Fig. 6. Schematic of the source separation estimation and usefulness of 
convex hull pixel selection for hyperspectral images. 
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